Diagnostic classification of cancers using DNA methylation of paracancerous tissues

The potential role of DNA methylation from paracancerous tissues in cancer diagnosis has not been explored until now. In this study, we built classification models using well-known machine learning models based on DNA methylation profiles of paracancerous tissues. We evaluated our methods on nine cancer datasets collected from The Cancer Genome Atlas (TCGA) and utilized fivefold cross-validation to assess the performance of models. Additionally, we performed gene ontology (GO) enrichment analysis on the basis of the significant CpG sites selected by feature importance scores of XGBoost model, aiming to identify biological pathways involved in cancer progression. We also exploited the XGBoost algorithm to classify cancer types using DNA methylation profiles of paracancerous tissues in external validation datasets. Comparative experiments suggested that XGBoost achieved better predictive performance than the other four machine learning methods in predicting cancer stage. GO enrichment analysis revealed key pathways involved, highlighting the importance of paracancerous tissues in cancer progression. Furthermore, XGBoost model can accurately classify nine different cancers from TCGA, and the feature sets selected by XGBoost can also effectively predict seven cancer types on independent GEO datasets. This study provided new insights into cancer diagnosis from an epigenetic perspective and may facilitate the development of personalized diagnosis and treatment strategies.

The potential role of DNA methylation from paracancerous tissues in cancer diagnosis has not been explored until now. In this study, we built classification models using well-known machine learning models based on DNA methylation profiles of paracancerous tissues. We evaluated our methods on nine cancer datasets collected from The Cancer Genome Atlas (TCGA) and utilized fivefold cross-validation to assess the performance of models. Additionally, we performed gene ontology (GO) enrichment analysis on the basis of the significant CpG sites selected by feature importance scores of XGBoost model, aiming to identify biological pathways involved in cancer progression. We also exploited the XGBoost algorithm to classify cancer types using DNA methylation profiles of paracancerous tissues in external validation datasets. Comparative experiments suggested that XGBoost achieved better predictive performance than the other four machine learning methods in predicting cancer stage. GO enrichment analysis revealed key pathways involved, highlighting the importance of paracancerous tissues in cancer progression. Furthermore, XGBoost model can accurately classify nine different cancers from TCGA, and the feature sets selected by XGBoost can also effectively predict seven cancer types on independent GEO datasets. This study provided new insights into cancer diagnosis from an epigenetic perspective and may facilitate the development of personalized diagnosis and treatment strategies.
Cancer continues to be a leading cause of mortality worldwide 1 . On the basis of the GLOBOCAN 2020 estimates of cancer incidence and mortality reported by the International Agency for Research on Cancer, an estimated 19.3 million new cancer cases (18.1 million excluding nonmelanoma skin cancer) and almost 10.0 million cancer deaths (9.9 million excluding nonmelanoma skin cancer) occurred in 2020 2 . Prediction of cancer stage and type as the important applications in cancer diagnosis are crucial for planning appropriate treatments.
Previous studies related to cancer diagnosis mainly focused on molecular data of tumor tissues 3-6 . Broët et al. proposed a new statistic for identifying gene expression features that detected tumor progression 4 . Rahimi et al. developed the highly time-efficient benders decomposition algorithm for the forest formulation (BDForest) to solve the problem of finding the similarity between different cancers, which is beneficial in classifying the stage of tumors 7 . Some studies have investigated molecular data of human pan-cancer and identified key biomarkers for prognosis and diagnosis of pan-cancer [8][9][10] .
Paracancer is the place where tumor growth and metastasis start. Wang et al. reported that the complement and angiogenesis pathways correlated with cancer progression were activated in the paracancerous tissues 11 . This finding revealed that the changes in paracancerous tissues are crucial complements to the conventional analysis of tumor tissues. Clinically, paracancerous tissues are more accessible to obtain than tumor tissues. Therefore, cancer stage and type prediction based on molecular data of paracancerous tissues may provide valuable information on understanding tumor stage progression and contribute to developing new approaches for cancer diagnosis. To the best of our knowledge, there have been no systematic studies that utilize data of paracancerous tissues to classify cancer stage and type.
In this study, we utilized extreme gradient boosting (XGBoost) 12 to discriminate tumor stage based on DNA methylation profiles of paracancerous tissues. The proposed XGBoost approach obtained better predictive performance than the other four machine learning methods. Moreover, our model extracted significant features from genome-wide DNA methylation profiles. GO enrichment analysis provided evidence that DNA methylation biomarkers of paracancerous tissues were closely associated with the progression of tumor stage. Additionally, we employed the XGBoost algorithm to build a multiclass classifier, which can accurately identify nine different cancers on the basis of DNA methylation profiles of paracancerous tissues from TCGA. The feature sets selected by the XGBoost model have high accuracy in cancer type prediction on independent GEO datasets.

Results
Predictive performance comparison for cancer stage prediction. For cancer stage prediction, we evaluated the predictive accuracy of five classification methods by fivefold-cv on nine datasets in this study. The area under the curve (AUC) of the receiver operating characteristic curve (ROC), the area under the precisionrecall curve (AUPR), accuracy (ACC), matthews correlation coefficient (MCC), Precision and Recall for different models on nine datasets were calculated as shown in Table 1. The AUC, AUPR, ACC, MCC, Precision, Recall of these algorithms obtained in each fold of all datasets can be found in Supplementary Table S1-S6, respectively. To compare the XGBoost model more intuitively with other machine learning methods, we further plotted ROC curves of five classification models, which were shown in Fig. 1.
We compare the AUC values of the five classification algorithms on nine datasets. We observe that XGBoost obtains significantly better results than baseline algorithm NB on all nine datasets, and it improves the average AUC by 22.6%. SVM outperforms NB on eight out of nine datasets and it improves the average AUC by 8.8%. RF performs better than NB on seven out of nine datasets and it improves the average AUC by 8.8%. KNN is able to outperform NB on all nine datasets and it improves the average AUC by 14.4%. We notice that the performance of NB is superior to the baseline (i.e. AUC = 0.5) on five out of nine datasets, further indicating that the other four classification models yield higher performance.
In addition, we calculate the ACC values of the five different classification models on nine datasets. Our results show that XGBoost predominantly outperforms NB, KNN, RF, and maintains comparable performance compared to SVM. XGBoost is superior to other models on six out of nine datasets and its average ACC is 0.694 that is 2.5%, 3.1% and 0.4% higher than NB, KNN and RF, respectively. Furthermore, XGBoost is still more competitive than other methods in terms of AUPR. The average MCC of XGBoost is 0.128, which is 88.2%, 80.3%, 13.3% and 21.9% higher than SVM, RF, KNN and NB, respectively. The average Precision of XGBoost is 0.385, which is 48.1% and 27.5% higher than SVM and RF, respectively. The average Recall of XGBoost is 0.399, which is 30%, 18.8% and 9% higher than SVM, RF and KNN, respectively.
In summary, these results clearly indicate that XGBoost achieves better performance on nine datasets by assessing AUC. Moreover, XGBoost marks the best results in most evaluation metrics, suggesting that it is superior to other classification methods. Specifically, the maximum AUC scores for nine datasets are 0. Gene ontology enrichment analysis. The GO analysis can identify biological pathways for revealing the relation between tumor progression and the CpG sites derived from the XGBoost model using the DNA methylation data of paracancerous tissues. GO analysis was conducted with the missMethyl package based on the top 10% significant CpG sites of nine cancers. Then we respectively obtained 147, 154, 155, 153, 151, 150, 152, 150, 152 important GO terms for KIRC, BRCA, THCA, HNSC, KIRP, LUSC, LIHC, COAD and UCEC. All GO enrichment results of nine datasets can be found in Supplementary Table S7-S15, respectively. From GO analysis, many enriched GO terms for each cancer were determined. However, it was difficult to analyze them one by one, and these GO terms may be redundant. Therefore, we elected to cluster them into more representative terms using the Cytoscape plugin ClueGO. ClueGO network diagram was visualized based on the following basic parameters: kappa coefficient was set to 0.1, three categories of GO were used for ontology files, where each node and line represented a term and the correlation between the terms, respectively. Different node colors denoted the classification of terms according to the functions. ClueGO network diagrams of nine cancers were depicted in Fig. 2. The results from ClueGO enrichment clearly illustrated that the most significant GO terms were transcription by RNA polymerase II, transmembrane receptor protein serine/threonine kinase signaling pathway, neuron projection guidance, female sex differentiation, DNA-binding transcription repressor activity, forebrain development and cell junction assembly in KIRC, BRCA, THCA, HNSC, KIRP, LUSC, LIHC, COAD and UCEC. Among the most significant GO terms, we determined common GO terms related to tumor progression in nine cancers. Further analyses for these GO terms can be found in the discussion.
Classification of distinct cancers using XGBoost. Based on DNA methylation data of paracancerous tissues, we exploited XGBoost to construct a tumor specific classifier for cancer type prediction. Table 2 presented all prediction results for TCGA and GEO datasets, respectively. XGBoost model obtained an accuracy of 100% for KIRC, 100% for BRCA, 100% for THCA, 100% for HNSC, 100% for KIRP, 100% for LUSC, 100% for LIHC, 100% for COAD and 100% for UCEC, exhibiting an average accuracy of 100% across all nine cancer www.nature.com/scientificreports/

Discussion
In previous studies, researchers mostly paid attention to the tumor itself and sought for the initiation and progression indicators of cancer from tumor tissues. However, more and more studies have suggested that inflammatory microenvironment is closely correlated with tumorigenesis and development. When the tumor develops to a certain stage, the paracancerous tissues are in the state of ischemia and hypoxia, and the increase of autophagy level will also promote the chemotherapy resistance, recurrence and metastasis of cancer, suggesting a poor prognosis 13 . In this sense, certain molecular responses and activities in paracancerous tissues may be related to the characteristics and status of cancer, thus providing potentially useful information for cancer type and stage prediction.
It is well recognized that reliable tumor stage prediction is critical for determining therapeutic strategies. Moreover, tumor treatment is highly dependent on the correct identification of the tumor origins. Tumor growth and metastasis affect paracancerous tissues, thereby it is valuable to investigate the relationship between molecular data of paracancerous tissues and pathological tumor stage.
To explore the potential role of DNA methylation profiles from paracancerous tissues in predicting cancer stages, we exploited the XGBoost algorithm to construct a classification model for cancer progression. We evaluated the performance of five classification methods by fivefold-cv.The comparison results of five advanced machine learning methods on nine datasets showed that XGBoost outperformed other classification models by assessing AUC and performed the best in the majority of metrics.
Previous studies almost utilized DNA methylation of tumor tissues for cancer stage prediction. Ma et al. employed the XGBoost model on the basis of multi-omics data to distinguish early and late stage tumors 14 . For KIRC, it achieved ACC of 0.719 and AUC of 0.797 based on DNA methylation data, whereas our results indicated that the ACC and AUC scores of XGBoost were 0.675 and 0.780, respectively. Deng et al. utilized gene expression and DNA methylation to build three networks for identifying the KIRC stages 15 . The prediction accuracy of the network using DNA methylation profiles was 0.696, which was 3.1% higher than the accuracy of our www.nature.com/scientificreports/    www.nature.com/scientificreports/ We performed GO analysis for the CpG sites identified by XGBoost model. The results indicate that the enriched GO terms associated with tumor progression are neuron projection guidance, cell junction assembly, transmembrane receptor protein serine/threonine kinase signaling pathway, transcription by RNA polymerase II and DNA-binding transcription repressor activity. Neuron projection guidance and cell junction assembly mainly refer to cellular processes. It is consistent with the knowledge that the defective functioning of cell biological processes is considered to be associated with tumor progression 17 . As a signal transduction pathway contributing to the pathogenesis of cancer, transmembrane receptor protein serine/threonine kinase signaling pathway describes a series of molecular signals initiated by the binding of an extracellular ligand to a receptor on the surface of the target cell 18 . Another GO term, transcription by RNA polymerase II, as the endpoint of signal transduction pathways, is the basis of development and differentiation 19 . The pathogenic mechanisms leading to cancer frequently involve altered signal transduction pathways. Furthermore, there remains a specific GO term of molecular function (MF), DNA-binding transcription repressor activity, which represses or decreases the transcription of specific gene sets. Aberrant regulation patterns at transcriptional level is regarded as a cause of human diseases 20 . Overall, these enriched GO terms demonstrated that underlying regulatory processes may participate in tumorigenesis. Moreover, the biological interpretation of the enriched GO terms suggested the relation between tumor progression and the significant CpG sites, which may serve as surrogate biomarkers for cancer diagnosis.
To our knowledge, it is the first research to establish a classification model for separating tumor stages on the basis of DNA methylation data of paracancerous tissues. Our study gave a systematic assessment of the performances of several machine learning algorithms for discriminating tumor stage and revealed the significance of paracancerous tissues for cancer progression. We also used XGBoost to construct the tumor specific multiclass classifier which can predict cancer type with high accuracy based on DNA methylation of paracancerous tissues. Furthermore, the utility of our model was emphasized by identification of the key CpG sites and GO terms associated with oncogenesis and tumor progression. Altogether, the investigation of DNA methylation profiles from paracancerous tissues may be helpful for understanding cancer progression and discovering new biomarkers. Our findings suggested that paracancerous tissues could be used as surrogate tissues for cancer stage prediction when tumor tissues were quite challenging to obtain.
Nevertheless, our study still has some limitations. First, the main limitation is the small sample size of nine cancer types used in the study. The inherent problem of the small sample size resulted in imprecision of prediction models. Second, we observed that the patients at different stages of KIRC and THCA can be well distinguished, whereas different stages of BRCA patients can't be clearly separated. The reasons for this are probably tumor heterogeneity and differences in tumor type. Because of tumor heterogeneity, there are obvious individual differences among BRCA patients 21 . Due to the imbalanced sample ratio of early stage and late stage being about 5:1, our model can not achieve better performance on LUSC. In addition, the reason for the low AUC of KIRP on our model may be due to a small sample size. Third, the races of patients include Asian, Black or African-American, White and not available, and White group accounts for the majority. Considering racial differences in cancers, our findings may not be suitable for paracancerous tissue data collected from other races.
In conclusion, our study suggested the potential role of paracancerous tissues in cancer diagnosis. One of our future efforts is to examine the possibility of other molecular data of paracancerous tissues in predicting the stage of tumors. The further application of our findings will contribute to understanding tumor progression and ultimately improving tumor treatment.

Methods
Data collection and pre-processing. We obtained DNA methylation profiles (HumanMethylation450; Level 3) and the corresponding clinical data of several cancers from The Cancer Genome Atlas (TCGA)database. Only paracancerous tissue samples were taken into account in the study. Paracancerous tissues in TCGA were For DNA methylation, we excluded the CpG sites with missing values in more than 20% of samples, and then imputed the remaining missing values using "na.roughfix" function in the "randomForest" package 24 . Table 3 presents the number of early stage samples, late stage samples and DNA methylation profiles for nine datasets.
We also collected clear cell renal cell carcinoma (GSE61441 25 ), breast carcinoma (GSE69914 26 ), thyroid carcinoma (GSE86961 27 ), head and neck squamous cell carcinoma (GSE75537 28 ), lung squamous cell carcinoma (GSE94785 29 ), liver hepatocellular carcinoma (GSE54503 30 ) and colon adenocarcinoma(GSE42752 31 ) from Gene Expression Omnibus (GEO). Each dataset contained paracancerous tissues and tumor tissues. The GEO datasets were only utilized for cancer type prediction due to the lack of pathological stage information. Therefore, we selected all paracancerous tissues of the seven datasets for independent validation. Table 4 shows the number of patients and DNA methylation profiles for the GEO datasets.
The samples from the public databases have obtained ethical approval. And all methods were conducted in accordance with the relevant guidelines and regulations.
Classification for cancer stage. We combined the samples annotated with stages I and II as early stage, and the samples annotated with stages III and IV as late stage. We used Python (version 3.7.3) and R (version   www.nature.com/scientificreports/ 4.0.4) for data pre-processing, classification and feature selection 32 . Figure 3 describes the framework developed for cancer stage prediction.
Extreme gradient boosting. Chen et al. developed a highly efficient and flexible gradient boosting algorithm called Extreme gradient boosting (XGBoost) 12 . It utilizes a more precise objective function and regularization term that improves its generalization ability 33 . This algorithm has been widely applied to many fields and shows its advantages in classification and regression studies [34][35][36][37] . Given a dataset D = x i , y i , here x i denotes CpG site of DNA methylation, y i ∈ {0, 1} is the class label of sample. Assuming that K is the number of trees, F represents the basic tree model. y i is denoted as the prediction of x i and given by Eq. (1): where f k (x i ) represents the predictive score of the k-th tree. Equation (2) denoted the objective function of XGBoost.
The loss function l denotes the mean square error between the prediction y i and the target y i . The term is utilized for smoothing the final learnt weights. The formula is shown as follows: where w denotes the score on each leaf, T denotes the number of leaves, γ and represents the degrees of regularization. The objective function at the t-th iteration can be described as Eq. (4).
where g i and h i refer to first and second order gradient statistics on the loss function.
(1) www.nature.com/scientificreports/ Other machine learning methods. For comprehensive analysis and comparison, we also employed other four well-known machine learning algorithms for building prediction models. Support vector machine (SVM) is a powerful supervised learning classifier 38,39 . As a kernel-based method, it aims to find optimal hyperplane that can perfectly distinguish different classes 40 . Random forest (RF) is a machine learning ensemble technique that constructs numerous decision trees based on different subsets of the data 41,42 . K-Nearest Neighbor (KNN) is a kind of simple classifier that has been extensively used for data classification 43 . Its performance is highly dependent on measuring the distance between the test samples and the training samples 44 . Naive Bayes (NB) is a probabilistic classifier that implements Bayesian techniques. The main characteristic of the classifier is that it's robust to noise and irrelevant attributes 45-47 . Model optimization. In this study, we employed fivefold cross validation (fivefold CV) for assessment of the classification models. The processes were that the dataset was randomly divided into five equal folds, and taken turns to use each fold to estimate the trained model, while the remaining four folds were used to train model. The relevant parameters for each model were optimized on the training set using threefold cross validation and grid search. The classification model was trained on the training set in combination with the optimal parameters. We obtained the performance metrics of the model by averaging all results of five test sets. For SVM, the parameters C and gamma were selected to optimize in the RBF kernel. For RF, the number of decision trees was adjusted. For KNN, we tuned hyperparameters of weights and the number of neighbors. For XGBoost, the configuration of parameters was a daunting task due to its many parameters. The optimized parameters included: 'learning_rate' , 'colsample_bytree' , 'subsample' , 'gamma' , 'min_child_weight' , 'max_depth' , 'reg_lambda' , 'reg_alpha' . Grid search made an exhaustive evaluation for various combinations of parameters and found the optimal set of parameters with the best performance.

Performance metrics.
To examine the performance of models, we utilized various evaluation metrics, commonly used to measure the classifier performances. The selected evaluation metrics include the area under the ROC curve (AUC), the area under the precision-recall curve (AUPR), accuracy (ACC), matthews correlation coefficient (MCC), Precision and Recall. FP, FN, TP and TN respectively indicate false positive, false negative, true positive and true negative predictions.
• AUC is applied to reflect the overall classification performance of the classifier by setting the discrimination threshold for comparing with the predicted probability from the classifier. • AUPR considers the recall and precision over different thresholds.
• Accuracy is metric of model robustness and represents the percentage of correct classifications by the classifier on the test set.
• MCC is commonly considered as a balanced indicator that can be utilized even though the classes are heavily imbalanced.
• Precision shows the ratio of correctly predicted positive samples accounts for the total number of the predicted positive samples.
• Recall shows the ratio of correctly predicted positive samples accounts for the total number of real positive samples.
Feature selection with XGBoost. We utilized the XGBoost algorithm to identify the key CpG sites that differentiated early-and late-stage cancers, and these features gave insight into the biological mechanisms of cancer formation and progression. For the XGBoost algorithm, the importance score of feature can be obtained on the basis of its participation in making key decisions with boosted decision trees 12 . All input features are ranked in descending order based on their importance scores. A higher score represents that the feature is more important. We selected the top 10% CpG sites as the significant feature sets to further explore the relationship between the feature sets and cancer stage.
Gene ontology enrichment analysis of the significant CpG sites. To explain the underlying biological mechanisms of above-mentioned important CpG sites identified by XGBoost, Gene Ontology (GO) enrichment analysis was conducted using gometh function in missMethyl package, taking into account the number of CpG sites per gene 48 . The GO terms with FDR < 0.05 were considered to be significant. Subsequently, Cytoscape www.nature.com/scientificreports/ (version 3.6.0) plugin ClueGO (version 2.5.7) (http:// apps. cytos cape. org/ apps/ cluego) was utilized to cluster GO terms and showed the distribution of the clusters over the GO terms, where the kappa statistic was set to greater than or equal to 0.1 49 . By using ClueGO, the redundant GO terms were reduced and the more representative terms were preserved in our study.

Construction of tumor specific multiclass classifier. Constructing a tumor specific classifier to
identify cancer type may be valuable in the common case where the primary origin of the tumor is unknown, since determining cancer type is critical to guide more detailed diagnosis and therapy. To this aim, we built a CpG-based tumor specific model using the XGBoost algorithm that can accurately classify cancer type. We first matched nine cancer datasets (KIRC, BRCA, THCA, HNSC, KIRP, LUSC, LIHC, COAD and UCEC) from TCGA with seven corresponding cancer datasets from GEO and retained a total of 208,745 common CpG sites for every dataset. We further merged DNA methylation profiles of different cancer types from TCGA and GEO, respectively. The sample sizes of TCGA and GEO datasets are 559 and 274, respectively. XGBoost model was built using 80% of the TCGA datasets as the training set, with the remaining 20% used to assess model performance. In addition, we utilized independent GEO datasets to validate our tumor specific classifier developed using TCGA datasets, which included thirty-one CpG sites as the feature sets by XGBoost-based feature selection. For multi-classification, we used accuracy to evaluate the performance of XGBoost.